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ABSTRACT 

We present a family of self-consistent, spherical, lowered isothermal models, consisting of 
one or more mass components, with parameterized prescriptions for the energy truncation 
and for the amount of radially biased pressure anisotropy. The models are particularly suited 
to describe the phase-space density of stars in tidally limited, mass-segregated star clusters in 
all stages of their life-cycle. The models extend a family of isotropic, single-mass models by 
Gomez-Leyton and Velazquez, of which the well-known Woolley, King and Wilson (in the 
non-rotating and isotropic limit) models are members. We derive analytic expressions for the 
density and velocity dispersion components in terms of potential and radius, and introduce a 
fast model solver in python (limepy), that can be used for data fitting or for generating discrete 
samples. 

Key words: methods: analytical - methods: numerical - stars: kinematics and dynamics - 
globular clusters: general - open clusters and associations: general - galaxies: star clusters: 
general 


1 INTRODUCTION 

The evolution of globular clusters (GCs) is the result of an interplay 
between stellar astrophysics (stellar and binary evolution, stellar 
mergers, etc.), dynamical two-body relaxation and the interaction 
with the tidal field of their host galaxy (Heggie & Hut 2003). De¬ 
spite this plethora of physical processes at work on their respective 
time-scales, the instantaneous surface brightness profiles and kine¬ 
matics of GCs are well described by relatively simple distribution 
function (DF) based models that need very few assumptions (Gunn 
& Griffin 1979; Meylan & Heggie 1997; Zocchi, Berlin & Varri 
2012 ). 

The relative simplicity of GC properties is owing to the ab¬ 
sence of gas and non-baryonic dark matter and the collisional na¬ 
ture of their evolution, which drives them to tractable properties, 
such as spherical symmetry, isotropy and (quasi-)equipartition be¬ 
tween different mass species (e.g. Spitzer 1987). Because the re¬ 
laxation time-scale of GCs is much longer than their dynamical 
time, their instantaneous properties can be described by models that 
satisfy the collisionless Boltzmann equation (see e.g. Chapter 8 in 
Berlin 2014). 

Two-body interactions in GCs evolve the velocity distribu¬ 
tion of stars towards a Maxwell-Boltzmann distribution, at least 
in the core, where the relaxation time-scale is short. Models with 
isothermal cores are therefore a good choice for fitting properties 
of GCs. An obvious starting point for a discussion on model choice 
is, therefore, the isothermal model. This model has an infinite spa- 
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tial extent and infinite mass (Chandrasekhar 1939) and to make the 
model applicable to real star clusters, the assumption of the ideal¬ 
ized Maxwell-Boltzmann distribution of velocities needs to be re¬ 
laxed. This can be done by changing the model such that stars have 
a finite escape velocity. Woolley (1954) developed such a model 
by simply ‘lowering’ the (specific) energy £ by a constant. The 
DF, which describes the density in six-dimensional phase-space as 
a function of E, is then simply f(E) = A exp[-(£ - 0(r,))/j^], for 
E < (f>(rt), and /(£) = 0 for £ > Here j is a velocity scale, 
which in the isothermal model equals the one-dimensional velocity 
dispersion and £ is reduced by the specific potential at the trun¬ 
cation radius r^, This truncation in energy mimics the role 
of tides due to the host galaxy, which makes it easier for stars to 
escape by reducing the escape velocity. The resulting models are 
nearly isothermal in the core, and have a finite mass and extent. 

The DF of these models is discontinuous at £ = ^(n). Michie 
(1963) and King (1966) avoided this by subtracting a constant from 
the DF introduced by Woolley, which makes the models continu¬ 
ous at E = Compared to the Woolley models, the density of 
stars near the escape energy is reduced in these models (hereafter 
referred to as King models), and they display a more gentle trun¬ 
cation of their density profile. The resulting, more extended, low- 
density envelopes make these models resemble real GCs more (for 
an in depth discussion on the effect of the truncation on the density 
profiles see Hunter 1977). The spherical, non-rotating limit of the 
models introduced by Wilson (1975), hereafter called Wilson mod¬ 
els, are models that are continuous both in the DF and its deriva¬ 
tive. This is achieved by subtracting an additional term linear in £ 
from the DF. These models are yet more spatially extended than 
King models. For some GCs in Local Group galaxies, the Wilson 
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models provide a better description of the observed surface bright¬ 
ness profiles compared to the King models (McLaughlin & van der 
Marel 2005; Carballo-Bello et al. 2012 also show that models that 
are more extended than King models better describe the surface 
brightness profiles of some GCs). 

An additional outcome of the two-body relaxation process is 
that it drives the velocity distribution of the stars towards isotropy. 
Isotropic models, defined by a DF that only depends on E, are 
therefore a natural choice for clusters that are in late stages of their 
evolution, near dissolution. At early phases, however, the velocity 
distribution in the outer parts is expected to be radially anisotropic. 
This is, first, because the (incomplete) violent relaxation process 
that takes place during their formation results in a halo of radial 
orbits (Lynden-Bell 1967). Secondly, two-body ejections from the 
dense core populate the halo with radial orbits on a two-body relax¬ 
ation time-scale (Spitzer & Shapiro 1972). Michie (1963) proposed 
a separable DF, dependent on E and on the (specific) angular mo¬ 
mentum J to introduce radial anisotropy (hereafter referred to as 
Michie-King models). The DF of the Michie-King models is the 
product of the isotropic DF with an exponential term with a de¬ 
pendent argument. This is similar to Eddington’s method of includ¬ 
ing radial anisotropy in the isothermal model (Eddington 1915). As 
a result, the inner parts of the models remain approximately isother¬ 
mal and isotropic, which is appropriate to GCs because there the re¬ 
laxation time is short, and anisotropy becomes important at larger 
distances from the centre. Near the truncation radius the models 
become isotropic again as a result of the energy truncation. The 
latter property has a somewhat coincidental resemblance to GCs, 
because near the Jacobi radius the orbits of stars gain angular mo¬ 
mentum due to the interaction with the (tri-axial) tidal potential (Oh 
& Lin 1992), therewith suppressing the amount of radial anisotropy 
near the truncation energy. A review of the effect of anisotropy on 
model properties can be found in Binney (1982). 

In real GCs, which contain multiple mass components, the re¬ 
laxation process drives the systems towards equipartition, result¬ 
ing in the heavier components being more centrally concentrated, 
a state which is often referred to as mass segregated. King models 
with different mass species were first introduced by Da Costa & 
Ereeman (1976) and have since been applied to take into account 
the effects of mass segregation in mass-modelling efforts of Galac¬ 
tic GCs (e.g. M3: Gunn & Griffin 1979, Omega Cen: Meylan 1987 
and larger samples of GCs: Pryor & Meylan 1993; Sollima et al. 
2012). Mass segregation is important for almost all of the Galac¬ 
tic GCs, given their short relaxation time-scales, relative to their 
ages (Henon 1961; Gieles, Heggie & Zhao 2011). Approximating 
multimass systems by single-mass models can lead to severe biases 
in the inferred properties of GCs (Shanahan & Gieles 2015; Sol¬ 
lima et al. 2015) and it is, therefore, desirable to have the ability to 
include multiple mass components in a dynamical model of a GC. 

It is our aim to develop a family of models that capture 
the general behaviour of collisional systems discussed above, and 
whose properties can be varied by parameters that can be con¬ 
strained by observational data. Davoust (1977) showed that the 
expressions for the DE of the isotropic Woolley, King and Wil¬ 
son models can be generalized by a DE in which the exponential 
function of E is reduced by the leading orders of its series expan¬ 
sion. This approach was further generalized by Gomez-Leyton & 
Velazquez (2014, hereafter GV14), who showed that solutions in 
between these models can be obtained (these models are briefly 
reviewed in Section 2.1.3). In this paper we extend the models 
of GV14 to allow for the presence of (radially biased) pressure 
anisotropy and multiple mass components. We present an efficient 


Poisson solver in python to facilitate the use of these models in fit¬ 
ting observational data, and in drawing samples from the models, 
which can be used as initial conditions for numerical simulations. 

The paper is organized as follows: in Section 2, we define 
the models and in Section 3, we illustrate their main properties. 
In Section 4, we present the code limepy' and our conclusions and 
a discussion are presented in Section 5. Supporting material can be 
found in the appendices. 


2 MODEL DEFINITION AND SCALING 

2.1 Single-mass models 

2.1.1 Distribution function (DF) 


The DF of the single-mass family of models is 


f(E, J^) - A exp 




E - 


( 1 ) 


for E < <^(rt), and 0 for £ > 4>(ri). The DF depends on two integrals 
of motion: the specific energy E = 12 + fir), with v the velocity 

and (j){r) the specific potential at distance r from the centre, and the 
norm of the specific angular momentum vector J = |rxv| = rv sin ??, 
where r and v are the position vector and velocity vector, respec¬ 
tively, and ?? is the angle between them. The energy E is lowered 
by the potential at the truncation radius f(ri). 

In equation (1) we introduced the function 


Ey(a, x) 


exp(x) a = 0 

exp(x)P(a, x) a > 0, 


( 2 ) 


where P(a,x) = y{a,x)IT(a) is the regularized lower incomplete 
gamma function (see Appendix D1 for the definition of this func¬ 
tion and its properties). Combining the exponential and the incom¬ 
plete gamma function into a single function Ey{a, x) has advantages 
in deriving the model properties (see GV14 and Appendix D2 for 
details on the behaviour of this function). 

A model is specified by three parameters: the central poten¬ 
tial, which is a required boundary condition for solving Poisson’s 
equation and defines how concentrated the model is; the anisotropy 
radius Ta, which determines the amount of anisotropy present in the 
system (for increasing ra the models are more isotropic); the trunca¬ 
tion parameter g, which controls the sharpness of the tmncation of 
the model (this parameter is called y in GV14). The physical units 
of a model are defined by two scales: the velocity scale s, and the 
normalization constant A, which sets the phase-space density and 
therewith the total mass M. For more information regarding scales 
and parameters of the models we refer the reader to Section 2.1.2. 

The isotropic models (ra — > oo) and their properties are dis¬ 
cussed in detail in GV14. For these models, and integer values of 
g, three well-known families of models are recovered: when g = 0 
we retrieve the Woolley (1954) models, for g = 1 we recover the 
King models (Michie 1963; King 1966), and for g = 2 we find the 
(isotropic, non-rotating) Wilson models (Wilson 1975)^. In prac¬ 
tice, the models defined by equation (1) are radially anisotropic for 
G ^ t\, because of the A dependence in the first exponential. When 


' LiMEPY is available from https://github. com/mgieles/limepy. 

^ The Woolley, King and Wilson DFs follow straightforwardly from equa¬ 
tions (1) and (2), because Ey(0,x) = exp(x), P(l,x) = 1 - exp(-x), such 
that Ey{l,x) = exp(x) - 1 and P(2, x) = I - exp(-x) - texp(-x), such that 
Ey{2, x) = exp(x) - I - X. 
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g = 1, the DF is the Michie-King model (Michie 1963), which is 
often used to fit GC data (e.g. Meylan 1987; Sollima et al. 2012). 

The potential is found by solving Poisson’s equation. For 
the self-consistent problem we consider here, the potential is com¬ 
pletely determined by the density p associated with the DF. This 
problem is non-linear, because the DF depends on the potential. 
Since the models defined by equation (1) are spherically symmet¬ 
ric, Poisson’s equation is 


dr \ dr / 


= 47rGp, 


(3) 


where the density is obtained by means of an integration of the DF 
over all velocities 

p^J d^yf(E,jf (4) 

In Sections 2.1.3 and 2.1.4, we derive analytic expressions forp as 
a function of (p and r. Note that only in the anisotropic case the de¬ 
pendence on the radial coordinate r is both implicit (through p, as 
in the isotropic case), and explicit, i.e. p{p, r). Having analytic ex¬ 
pressions for pip, r), avoids the need of solving a double integral at 
each radial step, making it significantly faster to obtain the solution 
to Poisson’s equation. In the next section we introduce a convenient 
set of units to solve the model. 


2.1.2 Scaling and units 


To solve Poisson’s equation, we use a dimensionless (positive) en¬ 
ergy E = p — k, with dimensionless potential p = ip(rt)-p)ls^, and 
k = v^/(2i^). As in King (1966), we consider the dimensionless 
density by normalizing p to its central value, i.e. p = p/po. In this 
way, Poisson’s equation in dimensionless form reads 


'A dr \ dr / 


-9p. 


(5) 


The dimensionless radius is now defined by the other scales: f = 
r/rs, with r^ = 9s^HAnGpn). This radial scale was introduced in 
King (1966) and is often referred to as the King radius. The factor 
of 9 was introduced to give r^ the meaning of a core radius, be¬ 
cause for models with moderately high central concentration, the 
projected density at rs is about one half of its central value. 

The Poisson equation is solved by assuming the boundary con¬ 
ditions at f = 0: 0 = 00 and d0/dr = 0. As mentioned in Sec¬ 
tion 2.1.1, the central potential 0o is one of the parameters that de¬ 
fine the modeF. 


2.1.3 Isotropic models 

We first briefly review the isotropic version of these models, as in¬ 
troduced by GV14. Many quantities can be derived from the DF. 
The density p is found by integrating the DF over all velocities 
(equation 4) and the pressure is found by taking the second veloc¬ 
ity moment of the DF 

p = (6) 

pO-2 = (7) 


^ This parameter is called Wq in King (1966). 

^ By considering the first velocity moment of the DF we find the mean 
velocity: for these models, this quantity vanishes eveiywhere. We also note 
that expressions for higher order moments of the velocity distribution can 
be derived, but these are beyond the scope of this paper. 


Here = 3cr^j is the mean-square velocity, cr[^ is the one¬ 
dimensional velocity dispersion and we introduce a dimensionless 
density integral and a dimensionless pressure integral {X^^~) 

JP ^ A r dkk'/^E.ig, p-k) = Eyig + |, 0), (8) 

ytt Jo 

I""' = ^ r dfc P^^E,(g, p-k) = 2,Eyig+\,f). (9) 

ytt Jo 

The results of these integrations follow straightforwardly from the 
convolution formula of the Eyia,x) function (equation Dll). An 
alternative derivation by means of fractional calculus is presented 
in Appendix B. The dimensionless density that appears in Pois¬ 
son’s equation is therefore p = 2 ?/Xq, where is the result of 
equation (8) evaluated at 0 = po. The dimensionless mean-square 
velocity is found from - a- j s^ = jlf. 


2.1.4 Anisotropic models 

Here we present the relevant quantities for the anisotropic case. 
The details of the derivations can be found in Appendix A, and the 
derivations by means of fractional calculus can be found in Ap¬ 
pendix B. To solve the anisotropic models, we introduce t = cost?, 
such that we can write the integral over the angles as 4;r f df. We 
further introduce p = rjr^ such that the density integral becomes 


die J' dt exp [kp^(t^ “ ^^^~Ey{g, p — k) 


s/n 

2 

fn 


f 




( 10 ) 


Here E(x) is Dawson’s integral and we refer to Appendix D3 for 
some properties of this function. To first order, E(x) oc x, and we 
thus find that for large i.e. small p, equation (10) converges to 
the integral of the isotropic model (equation 8). The solution of the 
integration gives as a function of p and p 


Eyig+fP) p^ 0^"ilFl(l,g+f,-#") 

l+p2 '^l-HpZ F(g-H|) 


Here iEi(a,b,x) is the confluent hypergeometric function whose 
properties are given in Appendix D4. For small p, the second term 
on the right-hand-side goes to zero and the solution converges to 
the isotropic result of equation (8). This expression for the density 
integral allows for fast computations of the right-hand-side of Pois¬ 
son’s equation and facilitates efficient solving of the anisotropic 
models. 

For the anisotropic models, we need to calculate both the ra¬ 
dial and the tangentiaF components of the pressure tensor, as well 
as the total pressure. The radial and tangential component of the 
velocity vector are defined as = v cos S and Vt = v sin 0 and for 


^ The tangential velocity comprises the two components vj = Vg+\J, where 
Vg = I'sinSsim/j and = I'sindcos^o. The corresponding components of 
the velocity dispersion tensor are equal to each other, and each of them 
accounts for half of the tangential component: cr^ = 2a-g = 2cr^. 
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the corresponding integrals we find 

= - I d^ I df exp[p^^(t^-l)l 0 - fc), (12) 

-VWJo Jo ^ ' 

= - I d^ I df exp[p^^(?^-l)l (1-- fc), 

V^Jo Jo ^ ' 

(13) 

= - I d^ I df exp[p^^(t^-l)l ^ - Q. (14) 

■VwJo Jo ' 


By carrying out these integrals as described in Appendices A2 and 
B2, we obtain 


XPd 


XPd 


X 


Ey{g + f ,^) 

1 + p2 1 + p2 

Ey(g+\,^) 2 

1 + (l+p2) 

lFi(l,g+ 

1 + 


^«^i|Fi(i,g+i-#^) 
r(g+i) 

2p^ ^^*2 

+ i+p2r(g+i) 

+ ifi(2,g + 1,-^p'^) , 


jpa- 


X 


£^r(.? + f,^) (3 + p^) p^ 

l+p2 (1 +p2) 1 +p2 P(g+ 7^ 

^J^iFi(l,g+3-0p2) + 2,Fi(2,g+3-0p2) 
1 + p2 ^ ^ 


(15) 


(16) 


(17) 


2 

Note that the expression for 2?°"' resembles the expression for 

of equation (11), in the sense that the functional form is the same, 

but all arguments and the power index that include g are increased 

2 

by 1. We already saw a similar resemblance between and 1'’“' 
in the isotropic case (equations 8 and 9, respectively). 

With these expressions for the density and pressure integrals, 
we defined most of the properties of these models that are of direct 
relevance for comparison to data. In Section 4.2, we discuss how 
the projected quantities can be derived. 


2.1.5 Limits 


In this section we consider some limits of the models. In the core, 
where p is small (r <K f.p, the model is isotropic. This is because 
the second terms in equations (11), (15), (16) and (17) vanish due 
to the multiplication by p^. 

Near the truncation radius the models behave like polytropes 
and are, therefore, also isotropic, because 

limiF|(l,a+l,-p^0) = 1, (18) 

^->0 

lim£^(fl,0)= —(19) 

^—>0 v{ci +1) 


and the p dependence disappears. In this regime, we find 


limp = 

^^0 


pS+3/2 

r(g + 5/2)’ 


lim pd"^ = 3 
^->0 


0S+5/2 

r(g + 7/2)’ 


limpirj = - limpiT^ = - limp(T^. 

0^0 3 0^0 2 0^0 


( 20 ) 

( 21 ) 

( 22 ) 


This suppression of the velocity anisotropy near the truncation ra¬ 
dius results naturally from the mathematical definition of the trun¬ 
cation, and is appropriate for tidally huncated systems (Oh & Lin 


1992). In A-body models a tangentially biased anisotropy is ob¬ 
served near r, (Sollima et al. 2015), which cannot be reproduced 
by the models presented here. However, it is likely that most of the 
stars with tangentially biased velocities are above the escape en¬ 
ergy, so-called potential escapers and these are not considered by 
these models, nor any other model we are aware off. 

Models with 0o 0 are close to pure polytropes over their 
entire radial range. In this regime, and for g = 7/2 (i.e. a polytropic 
index n = 5, equation 20), we recover the Plummer (1911) model, 
which is infinite in extent (p oc at large radii), but finite in mass. 
Polytropes with n > 5 (i.e. g > 7/2) are infinite in extent and will 
not be considered here. For g < 7/2 models can have a finite rt 
depending on both and ra (see GV14 and Section 3). 

In the cores of models with po » 0 the DF approaches the 
isothermal sphere, because 

lim Ey(a, ij)) = exp(^). (23) 

^—*oo 

Models with g —» oo also approach the isothermal sphere. To con¬ 
clude, these models approach the isothermal sphere in the limit of 
00 —> oo, independent of g, but also in the limit of g —> oo, inde¬ 
pendent of 00 - 


2.2 Multimass models 


It is possible to consider models with multiple mass components, 
by considering the DF as the sum of DFs of the form of equa¬ 
tion (1), each of which describes a different mass component with a 
mass-dependent velocity scale parameter. The first to do this were 
Da Costa & Freeman (1976), who calculated multimass King mod¬ 
els. For a multimass model with Acomp mass components, 2Acomp+2 
parameters are required in addition to the ones introduced in Sec¬ 
tion 2.1.1 for single-mass models. These additional parameters are 
the values for the component masses mj, the amount of mass in 
each component Mj, 6 and t]. The latter two parameters set the 
mass dependence of the velocity scale i j and the anisotropy radius 
of each component Paj, for which we adopt power-law relations 

(24) 

h.j = hii]- (25) 

Here fij = m^/m is the dimensionless mass of component j and in 
is the central density weighted mean-mass 


. YjjttijPoj 

m = — -. 

LjPoj 


(26) 


Note that in the multimass models, the values of s and are the 
velocity scale and anisotropy radius corresponding to in. The defi¬ 
nitions of 6 and 77 are such that the anisotropy profiles are approxi¬ 
mately mass independent when S = rj (see equation 1). The typical 
values considered for these parameters are (5=1/2 and 77 = 0. 

We notice that in the limit of infinite 0o the velocity scale ^j 
approaches the one-dimensional velocity dispersion of mass com¬ 
ponent j, CTid,;, hence the traditional assumption for (5 = 1/2 implies 
equipartition = fhcr\^j = constant). However, it is important 
to keep in mind that for multimass models with typical and realistic 
values of 0 o, the velocity dispersion of each component in the cen¬ 
tre is smaller than Sj and, therefore, there is no equipartition (see 
Section 3.2.1 and Merritt 1981; Miocchi 2006). 

To solve a multimass model self-consistently, we compute the 
density for each mass component as in equation (4) and add all 
components on the right-hand-side of Poisson’s equation. The de¬ 
tailed procedure is described in Gunn & Griffin (1979), and here 
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we only briefly summarize the required steps. The dimensionless 
Poisson equation to solve is 


j 


(27) 


where aj is the ratio of the central density of the 7 -th mass compo¬ 
nent to the total central density, such that 




(28) 


Poj lP(p'^^3>0, 0) 


(29) 


By considering multiple mass components, we introduce an eigen¬ 
value problem in the solution of Poisson’s equation, because the 
values of poj that yield the desired Mj values are not known a pri¬ 
ori. Therefore, as a first step to solve the model, we assume that 
Uj = Mj! Yij Mj, and we obtain the solution by iteration (see Sec¬ 
tion 4 for details). 


2.3 Normalization and potential energy 

In solving the models we have chosen to define the dimensionless 
quantities in terms of the density scale po and the velocity scale i- 
(Section 2.1.2). In some cases it is useful to have an expression for 
the normalization constant A in the DF (equation 1), for example, 
when fitting models to discrete data. From equation ( 6 ) we find that 
A relates to the other scales as 


A = 


Po 


(30) 


For the multimass models there is a normalisation for each com¬ 
ponent, Aj. The relation with the mass scale = MjM is = 
rjpo = where we introduced M - f pd^r. 

The total dimensionless (positive) gravitational energy U of 
the model is calculated from integrating the potential (King 1966) 


U = 



GM^ 

2ft 


(31) 


The second term has to be added because ^ is a lowered potential. 
Note that this integration of ^ over mass is readily obtained from 
solving Poisson’s equation. 


3 MODEL PROPERTIES 


3.1 Single-mass models 


3.J.J Density and velocity dispersion profiles 


In Fig. 1, we show the density profiles, the velocity dispersion 
profiles, and the anisotropy profiles for isotropic and anisotropic 
models with different values of the truncation parameter g. The 
anisotropy profile is computed from cr^ and cr^ as 




(32) 


In the case of isotropy y0 = O, 0<j8<l indicates radially biased 
anisotropy (with = 1 implying fully radial orbits) and for tangen¬ 
tially biased anisotropy <0. BecauseyS is a measure of anisotropy 
locally, we also quantify the total amount of anisotropy with 


K = 


2K, 


(33) 


$0 = 5 



logf 


Figure 1. Dimensionless density profile (top), velocity dispersion profile 
(middle) and anisotropy profile (bottom) for models with different trunca¬ 
tion parameters g (different colours). Isotropic models are shown with solid 
lines, anisotropic models with fa/fh = 1.5 with dashed lines. 


introduced by Polyachenko & Shukhman (1981). Here and Kt 
are the radial and tangential components of the kinetic energy, re¬ 
spectively. For isotropic models k = 1, and for radially biased 
anisotropic models /<■>!. Polyachenko & Shukhman (1981) found 
that for K > 1.7 ± 0.25 radial orbit instability occurs. We use this 
criterion to check the stability of the anisotropic models we calcu¬ 
late. 

In Fig. I we show anisotropic models characterized by fa/fi, = 
1.5. Because the (dimensionless) half-mass radius fj, is not known 
before solving the model, we find the value of f^ that gives the cor¬ 
rect ratio fa/fh iteratively. We see that all models are approximately 
isothermal in the centre. When increasing g, the models become 
more extended. Including radial anisotropy also results in a larger 
truncation radius. 

Note that, with this choice of fa/fh, the maximum value as¬ 
sumed by the anisotropy function for g = 0 (Woolley model) is 
about 0.4, while for g = 2 (Wilson model) it is possible to achieve 
j 8 - 1 in the outer parts of the model. This dependence of the maxi¬ 
mum value of yS on g does not imply that there are differences in the 
total amount of anisotropy: for all the anisotropic models shown in 
Fig. 1, indeed, we find k ^ 1.2. The ability to calculate models with 
more radial orbits (larger ft) without increasing the radial compo¬ 
nent of the total kinetic energy is important to keep in mind when 
considering other physical effects that can enhance or suppress the 
amount of radial orbits, such as the presence of a dark matter halo 
(Ibata et al. 2013) and the galactic tides (Oh & Lin 1992). In a 
forthcoming study, we quantify the presence of radial orbits in di¬ 
rect N-body models of tidally limited clusters (Zocchi et al. 2016). 
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3.1.2 DF, density of states and differential energy distribution 

In the top panels of Fig. 2, we show the DF as a function of E, for 
isotropic models, with different values of g and In the middle 
panels we show the density of states g(£), which is the phase-space 
volume per unit of energy (see equation C2 for a definition). The 
bottom panels display the differential energy distribution dM/d£, 
which is the amount of mass per unit energy. For the isotropic mod¬ 
els it is simply the product of f{E) and g(£’) (equation Cl). Details 
on how this was derived for the models presented here, and on the 
procedure for anisotropic models, are given in Appendix C. A gen¬ 
eral discussion on the differential energy distribution can be found 
in chapter 4 of Binney & Tremaine (1987). 

In the first and third columns (linear x-scale), we recognize 
the exponential behaviour of f{E) for the g = 0 model, and the 
exponential behaviour at high E for g > 0 models. From the second 
and fourth column, we see that at low E, the DF scales as f{E) oc 
£*, which corresponds to the regime where the models behave as 
polytropes. From Fig. 2 it is also evident that when E ^ fo, the 
model behaviour is independent of g. 

From the differential energy distribution, we see that only for 
g = 0 there is a non-zero mass at £ = 0. For models with g > 0, the 
truncation is such that f(E = 0) = dMjdE = 0. These models 
give rise to more realistic looking density profiles, but in real GCs 
the number of particles with the escape energy is not zero (Baum- 
gardt 2001), because of the gradual scattering of particles over the 
critical energy for escape by two-body relaxation, and because of 
the finite time for stars to escape from the Jacobi surface imposed 
by the galactic tidal field (Fukushige & Heggie 2000). 


3.1.3 Einite and infinite models 

As discussed in Section 2.1.5, there are no models with finite ex¬ 
tent if g > 3.5. GV14 showed that the maximum value gn,ax to get 
models with a finite extent depends on 0o, and g^ax = 3.5 holds in 
the limit of fa —» 0. GV14 show that all their isotropic models are 
finite for g < 2.1. 

We note that there is a class of isotropic models that are finite 
in extent, but are not relevant to star clusters, and that are not dis¬ 
cussed in GV14. This is illustrated in Fig. 3, where we show density 
profiles for models with different 0o and g = 2.75. The model with 
fo = 3 converges to a finite r, and has a density profile comparable 
in shape to the ones shown in Fig. 1. The model with ^g = 9 is in¬ 
finite in extent, and only plotted up to log f = 10. The models with 
0g = 5 and fo = 1 are finite, but show a sharp upturn in the density 
profile at large radii, which causes them to have a lot of mass in the 
envelope, but little energy, which makes these models inapplicable 
to real stellar systems. Their extreme density contrast between the 
core and the extended halo makes these models perhaps applicable 
to red giant stars (see the density profiles for red giants in Passy 
et al. 2012). To quantify the boundary between models with, and 
without the core-halo structure, we compute the ratio of the dimen¬ 
sionless virial radius = -GM^l(2U) over fh for a grid of models 
with 0 < ^g < 20 and 0 < g < 3.5, and we show the result as 
contours in Fig. 4. We find that for a given g(0o), when increasing 
0 g(g), the change in fv/fh is large and abrupt once the models de¬ 
velop the core-halo structure. We identify the value of fv/fi, 0.64 
as the one separating the two classes of models. In the remaining 
discussion, we only consider models with fv > 0.64rh. 

When considering anisotropic models, we find that for each 
f and g, there is a minimum value of that can be used to ob¬ 
tain a model that has a finite extent. We note that models with in¬ 


finite extent can have a finite total mass, but because we envision 
an application of these models to tidally limited systems we do not 
consider them here. In Fig. 5, we show the minimum fa for which 
models are finite in extent. The lines show, as a function of 0g, and 
for ditferent g, the values of fa that are needed to get f = lO’. Note 
that this minimum for fa goes up approximately exponentially with 
0 g, and also increases with g. 

3.1.4 Entropy 

King (1966) suggested that in the process of core collapse, clus¬ 
ters evolve along a sequence of models with increasing central con¬ 
centration. He also noted that his models are probably not able to 
describe the late stages of core collapse, because for large central 
concentration the variation in energy due to a change in the central 
concentration occurs in the envelope, and not in the core. Further 
support for this idea comes from Lynden-Bell & Wood (1968), who 
showed that a maximum in entropy occurs at ^g =: 9 for both Wool- 
ley and King models at constant mass and energy. The entropy of a 
self-gravitating system is obtained from the DF as 

S = - J'd^rdAflnf . (34) 

Because two-body encounters continuously increase the total en¬ 
tropy of the system, we do not expect King models to be able to 
describe a system in the late stages of core collapse (i.e. 0g > 9). 
This was confirmed by Fokker-Planck models of isolated star clus¬ 
ters going into core collapse (Cohn 1980), for which the entropy 
increase follows that of King models with increasing central con¬ 
centration, up to a value of ^g == 9, but then it continues to rise 
during the gravothermal catastrophe. Cohn concluded that in this 
regime, the isotropic King models are not able to describe the en¬ 
tropy evolution in his simulations. 

In Fig. 6, we show the entropy 5, computed as in equa¬ 
tion (34), for the isotropic King models (black solid line), which 
shows a maximum at 0g == 9. We also show the entropy curves for 
different values of g, and for selected anisotropic models. All mod¬ 
els are scaled to the same M and total energy Etot, in the conven¬ 
tional Henon A-body units: G = M = -4£,ot = 1 (Henon 1971). 
For 0 < ^g < 1, the anisotropic models are similar to their corre¬ 
sponding isotropic models, and therefore they have similar entropy. 
From this plot it is apparent that evolution at constant mass and en¬ 
ergy, and with increasing entropy is possible beyond 0 > 9 if g is 
increased, and/or is decreased (i.e. including more anisotropy). A 
local maximum in entropy is seen near fg ^ 17. Similar oscillating 
behaviour of the entropy was found for isothermal models in a non¬ 
conducting sphere and we refer to Lynden-Bell & Wood (1968) and 
Padmanabhan (1989) for detailed discussions. A study of equilib¬ 
ria in lowered isothermal models of the Woolley and King-type can 
be found in Katz et al. (1978); for a discussion on the evolutionary 
sequence of quasi-equilibrium states in A-body systems we refer to 
Taruya & Sakagami (2005). It would be of interest to compare the 
models discussed here to the phase-space density of particles in an 
A-body system undergoing core collapse. 

In Fig. 7, we illustrate the dependence of the entropy on g 
and 0g for isotropic models. For a model with g = 1 and a low 
concentration, the entropy can be increased by moving to the right 
in this diagram, and near fo - 9 the entropy can be increased by 
increasing g. 

In Fig. 8, we show the dependence of entropy on anisotropy, 
expressed here in terms of r^/rh, for models with g = 0. We see 
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Figure 2. Distribution function (DF, /(£)) (top), the density of states q(E) (middle) and the normalized differential energy distribution {dM/diE)/M = 
f{E)g(E)/M (bottom) for isotropic models with different central potentials (ipo = 3 in the two left columns and = 7 in the two right columns) and different 
g (see the legend in the top left panel). 


that for constant ra/ri, > 1, the entropy can increase by increasing 
00 , up to about 00-9 (this was also found by Magliocchetti et al. 
1998, in a study of anisotropic Woolley, King and Wilson models). 
The entropy can be increased further by decreasing the anisotropy 
radius. A maximum is found near 0o - 9 and - r^^. 

3.2 multimass models 

Multimass models with Acomp mass bins require, in addition to 
the parameters of the single-mass models, + 2 parameters 

(Section 2.2). There is, therefore, a large variety of models that 
can be considered, and many properties that we can chose from 
to illustrate the behaviour of these models. We decide to focus on 
two properties that highlight important features of these multimass 
models in relation to mass segregation. In a follow-up study (Peuten 
et al., in preparation) we present a detailed comparison between the 
multimass models and A-body simulations of clusters with differ¬ 
ent mass functions. 

3.2.1 The role of 6 

In Fig. 9 we show the dimensionless central velocity dispersion 
of each mass component, aujo, as a function of its mass nij for 
isotropic, 20-component models with different 0o and g. The mass 
bins are logarithmically spaced between 0.1 Mq and 1 Mq (note that 


the units are not important, because the model behaviour depends 
only on the dimensionless values mjlm), and Mj oc mf, which 
corresponds to a power-law mass function dNIdirij oc mf -^ (i.e. a 
GC-like mass function). The mass segregation parameter was set to 
d = 1/2 (for the definition of 6, see equation 24). 

Despite the fact that mjf is constant for all mass bins, there is 
no equipartition between the different mass species, i.e. does 
not scale as for the different mass components. This is be¬ 
cause only in the limit of infinite central concentration 0o —» oo, 
= CTidjo, but for realistic values of 00. the ratio (Tijjo/ij < 1. Be¬ 
cause the central potential for the lower mass components is smaller 
than the global 0o that defines the model, the truncation in ener¬ 
gies reduces a up more for low-mass components (Merritt 1981; 
Miocchi 2006). This is illustrated by the 0o = 16 model in Fig. 9, 
for which a constant only holds for the most massive bins. 

Trenti & van der Marel (2013) recently observed very similar 
trends between crujo and mj in A-body models of GCs (see their 
fig. 1) as those shown in Fig. 9. They concluded that modelling 
techniques that assume equipartition, such as multimass Michie- 
King models, are ‘approximate at best’. We stress that multimass 
models that are widely used in literature, i.e. those with 6=1/2 (Da 
Costa & Freeman 1976; Gunn & Griffin 1979), are not in a state of 
equipartition, as is illustrated in Fig. 9 and has been stated previ¬ 
ously (Merritt 1981; Miocchi 2006). In fact, from a comparison of 
the model behaviour in Fig. 9 and the A-body models of Trenti and 
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Figure 3. Density profiles for isotropic models with truncation parameter 
g = 2.75. Models with = 3, 5, and 7 (blue, green, and red line, respec¬ 
tively) have a finite truncation radius, but only the model with = 3 is 
relevant when describing GCs; the model with = 9 (light blue line) is 
infinite in extent. 



$0 


Figure 4. Ratio of dimensionless virial radius to half-mass radius, ry/rh^ for 
models with different 0o and g. We consider models with ry/rh ^ 0-64 as 
relevant to describe star clusters. Models that have an infinite rt are plotted 
as Ty/rh = 0 (i.e. they coiTespond to the white region in the plot). 


van der Marel we conclude that the most commonly chosen flavour 
of multimass models (i.e. King models with (5=1 /2) do a good job 
in reproducing the degree of mass segregation in evolved stellar 
system (see also Sollima et al. 2015). 



Figure 5. Minimum ra for finite sized models, for different and g. 



Figure 6. Entropy curves for isotropic and anisotropic models with differ¬ 
ent truncation prescriptions (i.e. different values of g). All models are scaled 
to the same mass and energy. The anisotropic models are shown as dashed 
lines and for these models we used r-^ = r^. For g > 1.5 the anisotropic mod¬ 
els are not finite for all and the corresponding curves are therefore not 
plotted. This figure shows that the entropy can be increased by increasing 
g, and/or by decreasing ra. 


i.2.2 The role ofrj 

In Fig. 10, we illustrate the effect of the parameter rj that sets the 
anisotropy radius of the different mass components (for the defi¬ 
nition of rj see equation 25). We show the anisotropy profiles for 
three-component models with rrtj = [0.2,0.4,0.8], and the same 
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Figure 7. Entropy contours for isotropic models, all scaled to G = M = 
-4iitot = 1: with different and g. Contours of constant ft are shown as 
black lines. Moving to the right at constant g leads to an increase of entropy 
up to ^0 - 9 (Lynden-Bell & Wood 1968; Cohn 1980). The entropy can 
grow further by increasing g at constant - 9. The maximum entropy is 
found near - 9 and g ^ 2.2. 



Figure 8. Entropy for models, all scaled to G = M = -4iitot = 1 ■ with g = 0 
and different concentrations and different amounts of anisotropy, quantified 
here in terms of ra/rj,. Contours of constant ft are shown as black lines. A 
maximum in entropy is found at ^ 9 and ra ^ Hi. 


mass function as before (i.e. Mj oc mf), and for different values of 
Tj. All models have = 9, g = 1.5, d = 1/2 and fa = 20. 

In the multimass models used in the literature 7/ is implicitly 
assumed to be 0. From Fig. 10 we can see that this implies that the ft 
profile of the high-mass stars rises to larger values. It is tempting to 
interpret this as that massive stars are on more radial orbits. How¬ 
ever, the more massive stars are also more centrally concentrated, 
where the velocity distribution is more isotropic. To quantify the 



Figure 9. Dimensionless central velocity dispersion for each mass compo¬ 
nent of multimass models with a power-law mass function, and different 0o 
and g. The value of the mass segregation parameter is d = 1/2. Equipar- 
tition in energy is only reached for large values of for the model with 
00 = 16. The dash-dotted lines show the velocity dispersion each of the 
models would have in the case of equipartition. 


importance of this effect, we show in each panel the values of the 
parameter kj for each mass component (equation 33). From this we 
can see that in fact for the rj = 0 models the intermediate mass 
component is the most anisotropic. The relation between and kj 
depends on the mass function, 0o> and g and this is therefore not a 
general property of = 0 models. 

We note that for 77 = d = 1/2 the Pj profiles are nearly mass 
independent. Again, this does not mean that the kinetic energy in 
radial orbits relative to that in tangential orbits is constant, as can 
be seen from the values of kj. When considering a value of 7 / > 1/2 
we observe that the component for which Pj assumes the largest 
values is the least massive one. 


4 THE LIMEPY CODE 
4.1 General implementation 

We introduce a PYTHON-based code that solves the models and al¬ 
lows the user to compute some useful quantities from the DF. The 
code is called Lowered Isothermal Model Explorer in PYthon 
(limepy), and is available from: https: //github. com/mgieles/ 
limepy. 

One of the main features of the code is its flexibility: the user 
can easily solve isotropic or anisotropic models, and include one 
or more mass components. The type of model to calculate is deter¬ 
mined by the input parameters: 

(i) the dimensionless central potential 0 o; 

(ii) the truncation parameter g; 

(iii) the anisotropy radius f.^ (for anisotropic models); 

(iv) two arrays nij, Mj and 6 and 77 (for multimass models). 


MNRAS 454, 576-592 (2015) 



























10 


Mark Gieles and Alice Zocchi 



logr 


Figure 10. Anisotropy profiles for three-component models (nij = 
[0.2,0.4,0.8]) with dilferent values for the anisotropy parameter rj, that sets 
the anisotropy radius of the individual components as a function of their 
mass. All models have 0o = 9, g = 1.5, i5 = 1/2 and = 20. The values of 
Kj are shown for each component in the individual panels. 


By default, the model is solved in the dimensionless units de¬ 
scribed in Section 2.1.2. There we pointed out that the scales of 
the models are set by A and s, which correspond to a mass density 
(in six-dimensional phase space) and a velocity scale. These two 
scales, combined with the gravitational constant G then define the 
radial scale. To allow a user to scale a model to physical units, we 
decided to use the mass and radial scale as input, and from this the 
velocity scale is computed internally. The reason for this is that we 
foresee that an important application of the code is to recover the 
GC mass and radius from a comparison of the models to data. It 
is possible to scale the model to physical units by specifying M in 
Mq and a radial scale (either r, or Th) in pc. The resulting unit of 
velocity is then kms^*, with G = 0.004302pc(kms“')^/MQ. Al¬ 
ternative units, such as the Henon units G = = M = I (Henon 

1971), can be considered by redefining the scales. After solving the 
model, the values of all typical radii are available: the King radius 
rs, the half-mass radius Th, the truncation radius n, the anisotropy 
radius r^, and the virial radius r,. 

The code solves Poisson’s equation with the ‘dopriS’ integra¬ 
tor (Hairer, Nprsett & Wanner 1993), which is a Runge-Kutta in¬ 
tegrator with adaptive step-size to calculate fourth and fifth order 
accurate solutions. It is supplied by the scipy sub-package inte¬ 
grate. The relative and absolute accuracy parameters are chosen 
as a compromise between speed and accuracy and can be adjusted 
by the user. The basic version of the code allows us to obtain, as a 
result of this integration, only the potential as a function of radius. 
The full model solution contains, in addition to the potential, the 
density, the radial and tangential components of the velocity dis¬ 
persion, the global velocity dispersion profile, and the anisotropy 
profile (equation 32). It is possible to use the potential calculated 
in this way to compute the value of the DF as a function of input 


E (isotropic models), or E and J (anisotropic models), or positions 
and velocities. 

After solving a model, the code carries out a simple test to 
see whether it is in virial equilibrium: 2K - U = 0, where K is 
the dimensionless total kinetic energy (recall that U is defined to 
be positive, equation 31). For models that are infinite in extent, the 
solver stops at a large radius, the virial equilibrium assertion fails 
and the lack of convergence is flagged. 

For multimass models the central densities of the components 
need to be found by iteration (Section 3.2). Gunn & Griffin (1979) 
proposed a recipe in which the ratios of central densities over the 
total density, aj, are set equal to Mjj 2^ in the first iteration. Be¬ 
cause for d > 0 the more massive components are more centrally 
concentrated, the amount of mass in these components is under¬ 
estimated in the first iteration, while the mass in low-mass stars 
is overestimated. After each iteration, aj is multiplied by the ratio 
of MjIMj, where M'. is the array of masses obtained in the previ¬ 
ous step, and then normalized again. This is repeated until conver¬ 
gence. However, we found that for models with low and a wide 
mass spectrum the mass function does not always converge with 
this method. We found that multiplying ajhy instead, 

is more reliable and results in a similar number of iterations for 
models that do converge with the method proposed by Gunn and 
Griffin. 

Solutions are not numerically stable when considering large 
values of the arguments of the hypergeometric functions. To stabi¬ 
lize the calculations, we adopt the asymptotic behaviour of the hy¬ 
pergeometric function iFi(l,b,-x) and iFi(2,b,-x) for x > 700 
(see equations D24 and D25). For multimass models with a wide 
mass spectrum (e.g. when stellar mass black holes are considered 
in addition to the stellar mass function), the central potential of 
the massive component can be too large for the computation of 
f) (see equation 29) in the first iteration. We therefore use 
the approximation pj = — ^o)] if > 700. 

4.2 Model properties in projection 

In order to compare the models to observations of GCs, it is nec¬ 
essary to compute the model properties in projection. For a spher¬ 
ically symmetric system, it is straightforward to compute the pro¬ 
jected properties as a function of the projected radial coordinate 
R (for a more detailed discussion, see for example Binney & 
Tremaine 1987). 

The projected surface mass density is found from the intrinsic 
mass density as 

E(R) = 2 f‘dzp(r), (35) 

Jo 

where d = + d, and z is along the direction of the line-of- 

sight. The velocity dispersion along the line of sight is given by the 
following integral 

o'losW = ^ dzp(r)cr^(r), (36) 

where crj is the contribution of the velocity dispersion tensor to 
the z-direction. For isotropic models, cr^ = cr^/3. For anisotropic 
models, it is possible to calculate 

cr^(r) - aj cos^ ^ + crl sin^ (37) 

where sin^ = R/r. We recall that, for the anisotropic models con¬ 
sidered here, cr^ = cr^ = 0-^12. The component cr^ does not con¬ 
tribute to <T^, because it is always perpendicular to the line of sight. 
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We can use the anisotropy profile f} (see equation 32) to rewrite 
equation (37) as 


0-2 (r) = 0-2 


l-Pd)- 


(38) 


The quantity o-los( 7?) is useful when comparing the models to 
the velocity dispersion profiles that are calculated from radial (i.e. 
line-of-sight) velocities. Now that proper motions data are becom¬ 
ing available for an increasing number of GCs (Bellini et al. 2014), 
it is also interesting to compare the velocity dispersion components 
that can be measured on the plane of the sky with those calculated 
from the models. This comparison is particularly important because 
it is a direct way to detect the presence of anisotropy in the systems. 
We calculate, therefore, the radial and tangential projected compo¬ 
nents of the velocity dispersion as 


o-r(R) = ^ dzp(r)cTl(r), 
dzp(r)rr2(r), 
where cTj is given by 


(Tsd) = o-f 


1 -P(r) 


1 - 7?2 y 


(39) 

(40) 

(41) 


In the case of multimass models, the projected quantities in¬ 
troduced above are calculated separately for each mass component, 
by replacing every quantity in the equations above with the respec¬ 
tive yth profile. 


4.3 Generating discrete samples from the DF 

A separate sampling routine limepy. sample is provided that gener¬ 
ates discrete samples from the models. The routine takes a python 
object containing a model as input and the number of points N that 
need to be sampled. In the case of a multimass models the input 
N is ignored and computed from the total mass M and the pair 
nij, Mj. Radial positions are sampled by generating numbers be¬ 
tween 0 and 1 and interpolating the corresponding r values from 
the (normalized) cumulative mass profile(s). 

To obtain velocities, we first sample values of x, where x = 
p /2 _ (v 2/2)2/2. The probability density function (PDF) for x can 
be written as 

P(X) = 0(f) - x2/2), (42) 

where p = rlr,,. The function P(x) has a maximum at a = 0, and 
declines monotonically to 0 at x = 0(f)2/2. These properties make 
it easier to efficiently sample values for x from P(x), than sampling 
values of v from df{r,v). To make the rejection sampling more 
efficient, we adapt a supremum function F{x), which consists of 10 
segments between 11 values x, which are linearly spaced between 
0 and 0(f)2/2, and for each segment x, < x < x,+i, F,(x) = P(Xi). 
We then sample values from the function F{x), reject the points that 
are above P{x) and resample the rejected points until all points are 
accepted. Typically, a handful of iterations are needed. 

For anisotropic models we also need to sample angles 0. We 
do this by sampling values for t = cos 9. From the DF it follows 
that the PDF for t is 

P(f) = exp[p2fc(f2_i)]. (43) 

By integrating equation (43) we find that the cumulative DF is the 


imaginary error function. This function cannot be inverted analyti¬ 
cally, hence the values for t need to be found by numerically invert¬ 
ing this function, which can be done accurately with built-in scipy 
routines. 

When values for r, v, and are obtained, these are converted 
to Cartesian coordinates by generating three additional random an¬ 
gles. 


5 CONCLUSIONS AND DISCUSSION 

In this study we present a family of lowered isothermal models, 
with the ability to consider multiple mass components and a vari¬ 
able amount of radially biased pressure anisotropy. The models ex¬ 
tend the single-mass family of isotropic models recently developed 
by GV14. The new additions we propose here make the models 
ideally suited to be compared to data of resolved GCs. 

The models are characterized by an isothermal and isotropic 
core, and a polytropic halo. The shape of the halo is set by the 
truncation parameter g, that controls the sharpness of the energy 
truncation, i.e. the prescription of lowering the isothermal model. 
For integer values of g, several well-known isotropic models are 
recovered: for g = 0 we recover the Woolley (1954) models, for 
g = 1 the Michie (1963), or King (1966) models and g - 2 cor¬ 
responds to the non-rotating, isotropic Wilson (1975) models. The 
DF proposed by GV14, with the introduction of the continuous pa¬ 
rameter g to determine the truncation, allows us to consider models 
in between these models. The advantage of this prescription for the 
truncation is that it is now possible to control the sharpness of the 
truncation by means of a parameter. 

We present Lowered Isothermal Model Explorer in PYthon 
(limepy), a PYTHON-based code that solves the models, and com¬ 
putes observable quantities such as the density and velocity disper¬ 
sion profile in projection. In addition, the code can be used to draw 
random positions and velocities from the DF, which can be used to 
generate initial conditions for numerical simulations. 

It is interesting to discuss possible extensions of, and improve¬ 
ments to the models. One obvious pitfall is that the tidal field is not 
included in a self-consistent way. To quantify the effect of the omis¬ 
sion of the tidal field, we can consider the specific energy E at Xf 
In our models E(ri) = 0(rt) = -GM/rt, whilst inclusion of the tidal 
terms would give (for a cluster on a circular orbit, in a reference 
frame corotating with the galactic orbit) a specific Jacobi energy 
of £j(rt) = -(3/2)GM/rt. Therefore, the properties of stars near 
the critical energy for escape are described only approximately by 
these models, because in this energy range the galactic tidal poten¬ 
tial is of comparable importance as the cluster potential. Another 
simplification of the models is that the galactic tidal potential is tri- 
axial, whilst our models are spherical. Both of these points could be 
improved upon by including a galactic tidal potential in the solution 
of Poisson’s equation, following the methods described by Reg¬ 
gie & Ramamani (1995) or Berlin & Varri (2008), Varri & Berlin 
(2009). 

The models do not include a prescription for rotation, which 
can be an important factor to take into account when describing 
real GCs (e.g. Bellazzini et al. 2012). Self-consistent models with 
realistic rotation curves exist (Varri & Berlin 2012) and have been 
successful in describing the rotational properties of several Galactic 
GCs (Bianchini et al. 2013). It is feasible to include rotation in the 
models presented in this paper, for example, in the way it is done 
in the Wilson (1975) model, by multiplying the DF in equation (1) 
by a dependent exponential term. Including the rotation, and a 
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description of the galactic tidal field, would make the models more 
realistic and, therefore, a worthwhile exercise for future studies. 

Lastly, we note that our models could be useful in modelling 
nuclear star clusters. Despite the fact that these systems are not 
tidally truncated in the same way as clusters on an orbit around the 
galaxy centre, their profiles are well described by lowered isother¬ 
mal models (e.g. Georgiev & Boker 2014). For a general appli¬ 
cation to nuclear star clusters, it is desirable to include the effect 
of the presence of a black hole in the centre, which generates a 
point-mass potential. Miocchi (2010) provided a method to self- 
consistently solve King models with an external point-mass poten¬ 
tial: this recipe could be used to include the effect of a massive 
black hole in the models presented here, to make them more versa¬ 
tile in describing nuclear star clusters. 

The aim of this project was to introduce models that can be 
used to describe the phase-space density of stars in tidally limited, 
mass-segregated star clusters, in any stage of their life-cycle. At 
early stage, GCs are dense with respect to their tidal density (e.g. 
Alexander & Gieles 2013) and at the present day about half of the 
GCs is still much denser than their tidal density (Baumgardt et al. 
2010; Gieles et al. 2011). These GCs ought to have a population 
of stars with radial orbits in their envelopes, either as a left-over of 
the violent relaxation process during their formation (Lynden-Bell 
1967), and/or because of two-body ejections from the core (Spitzer 
& Shapiro 1972). In this phase we expect models with high val¬ 
ues of g, and small to describe GCs well. These models can thus 
describe GCs with large Jacobi radii, relative to rn. This applies to 
a large fraction of the Milky Way GC population, and these ob¬ 
jects are beyond the reach of King models (Baumgardt et al. 2010). 
In later stages of evolution, GCs will be more tidally limited, and 
isotropic, hence we expect g to reduce and r.^ to increase during 
the evolution (up to a value that practically corresponds to hav¬ 
ing isotropic models). Capturing these variations in GCs properties 
with continuous parameters has the advantage that these parameters 
can be inferred from data. This avoids the need of a comparison of 
goodness-of-fit parameters of different models. 

When only surface brightness data are available, it is challeng¬ 
ing to distinguish between models with different truncation flavours 
and pressure anisotropy, because their role has an impact mostly on 
the low-density outer parts, far from the centre of the cluster, where 
foreground stars and background stars are dominating. The addition 
of kinematical data of stars in the outer region of GCs greatly aids 
in discriminating between models, but this is challenging at present. 
Precise proper motions (< 1 km s“') can be obtained with the Hub¬ 
ble Space Telescope (HST; e.g. McLaughlin et al. 2006; Watkins 
et al. 2015), but the field of view of //ST limits observations to the 
central parts of Milky Way GCs. Radial velocity measurements of 
stars in the outer parts of GCs are expensive because of the contami¬ 
nation of non-member stars (Da Costa 2012). The upcoming data of 
the ESA-Gaifl mission will improve this situation: the availability 
of all-sky proper motions and photometry measurements will facil¬ 
itate membership selection, and for several nearby GC the proper 
motions will be of sufficient quality that they can be used for dy¬ 
namical modelling and to unveil the properties of the hidden low- 
energy stars (An et al. 2012; Pancino et al. 2013; Sollima et al. 
2015). The models presented in this paper allow for higher level of 
inference of physical properties of GCs from these upcoming data. 

In two forthcoming studies we will compare the family of 
models to a series of direct A-body simulations of the long term 
evolution of single-mass star clusters (Zocchi et al. 2016) and mul¬ 
timass clusters (Peuten et al., in preparation) evolving in a tidal 
field. 
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APPENDIX A: DERIVATIONS 

The DF introduced in equation (1) can be expressed as a function 
of the dimensionless quantities f k and p defined in Sections 2.1.2 
and 2.1.4 as 

/ = A exp sin^ 0 ) ^ - fc). (Al) 

We want to calculate, for these models, the density and velocity dis¬ 
persion components. We recall that these quantities can be obtained 


from the DF in the following way:® 


P = 

0-? = ^ J' d\fvf. 


(A2) 

(A3) 


where the subscript i denotes the i-th component of the veloc¬ 
ity vector. To carry out these integrals of the DF in the three- 
dimensional velocity volume we can use the dimensionless variable 
k and the variable t = cos 9 


d^v = dvdSdi^v^ sin0 = —dkdtAip 


(A4) 


In calculating the relevant quantities mentioned above, we en¬ 
counter the following integrals with respect to the variable t 


f df exp |^-fcp^(l - r“)j = 
•do 

J' dtd exp [-^p^(l - f^)j 


flp 

1 

2( flpf' 


(A5) 

(A6) 


where F{x) is the Dawson integral, whose properties are presented 
in Section D3. We use the above results to proceed and derive the 
density and velocity components. 


AI Density profile 

The density is calculated as 


f 


d\f 


A 

f/a 


— die J' dt exp - l)j (g, 0 - p) 
dk ^ ^ ' E,[g,cl,-k) 


Vtr 
: Ai'’, 


(A7) 


where we replaced F(3/2) by -0r/2, and we introduced A = 

I ,\ 3/2 

A127rj 1 and we solved the integral over t as shown in equa¬ 
tion (A5). The integral 2? can be solved by first doing an integra¬ 
tion by parts (by using the results in equations DIO and D15) and 
by then using the convolution formula of equation (Dll) and the 
recurrence relation of equation (D8) in the following way: 


IP =2 

W Jo ' ' 


(A8) 


^~^Jo 

= Ey(g+\,f)-p^— \ dkEy(g,^-k) -;- 


+ p 

yTT 


Jo 


dk 




r(g+i) 


- Ey 


g+lf)-p^T + p^TfX- 


(A9) 


® As noticed already in Section 2.1.3, equation (A3) holds because the 
mean velocity for the systems described by these models is zero every¬ 
where. 
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The integral Tp i can be calculated by substituting the Dawson func¬ 
tion for its series representation (see equation D14), by changing 
variable to y = kjij), by using the Beta function of equation (D12), 
and by recognizing the expression in equation (D16) 


2 

VIp) 

-0r Jo 
0s+iiEi( 

r(g+i) 

1..?+ 

p 

’) 

T 

i) 


(AlO) 


where iFi (a,b,x) is the confluent hypergeometric function (see 
Section D4). Therefore, we can finally write the density integral 
as 


Ey{g+l,?>) 0^+JiFi(l,g-l- 


(Alll 


A2 Velocity dispersion profiles 

The velocity dispersion profile can be computed in a similar way 
as the density, by using again the result in equation (A5) 




dAv^f 

2As^ 

yfn 

^ I - x-^kEy(g, ^-k) 

yn 


r’ . r r - n - - - 

— I dk dt e.x\)\-kp^{l - d)\jd^^Ey{g,^-k) 
^ Jo Jo ^ 

f 


(A12) 


2 

The integral TT”" can be solved with an integration by parts, then 
by using equation (Dll), and finally, after having used the recur¬ 
rence relation of equation (D8), by recognizing the presence of the 
integral found when calculating the density 





dk kEy (g, ^ - kj 




1,0 — kj 



^Ey(g+l,0) + 22:;^, + (Xp,2 -x^^). 


(A13) 


(A14) 


where = X^(g -l- l,p, 0). We then solve the integral Xf 2 in a 
similar way as we did forXfj, to get 



+2,Fi (2,g+l-p^)]. 


(A15) 


Therefore, we finally have 
As^ 


cdp 


(p2 -I- 1) 

r(g+i) 


Ey[g+h$)\j 


2 + 


+ P' 


[2,F, (2,g+^,-p2 




,Fi(l,g+|,-p2- 


2 + p^ 




(A16) 


The radial component of the velocity dispersion is given by 


rrp-/ 


d^v(v COS0)^/ 


I dk'iJI^Eyi^,0-k) 

VTT Jo ^ 


kp^ (yflpP 




— [Ey{g^\,0)-V], 


(A17) 


where we solved the integral by using equations (A6) and (Dll). 
We point out that we can express this quantity by means of the 
density integrals of the isotropic and the anisotropic case (see also 
equations 8 and 10 in Section 2). We finally obtain 

Ey{g+\,0) 0’^*^lFi{l,g+\,-p^0) 


2 r 7 

(T^P = As- 


(i+p2) (i+p2)r(g+|^ 


(A18) 


which, by using equations (D8) and (D19), can be rewritten as 


a-^p=As^ 




(1+P^) 


(l+p2)r(g+7) 


(A19) 


To calculate the tangential component of the velocity dispersion we 
solve the integral over t by expressing it as the difference between 
equation (A5) and equation (A6), and we carry out an integration 
by parts 


r;p=/ 


d3v(v sin0)V 

p d^ p/^Ey (g, 0-k) 
yn Jo ^ 


yfn 

2F I dkp 


1 F[dkp) 

+ 


yfkp kp^ yfkpP 


^[p^X^^-Ey{g+l,4:}+X^]. 


(A20) 


After recognizing the integrals we solved above, we can finally 
write 


cdp : 


A A 

(1 p2) 

r(g+i) 


[Ey[g+l$]j^ 


PJ 


iFi(l,g-l- l,-0p^) 


(l+p2) 


iFi (2,g-l- l,-0p^)]] ■ 


(A21) 


APPENDIX B: DERIVATIONS USING ERACTIONAL 
CALCULUS 

Fractional calculus is a branch of mathematics that considers real 
numbers for the orders of derivatives and integration. Because the 
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integrals that need to be solved contain terms like and we 
can use semi-derivatives and semi-integrals and integration by parts 
to solve them. By following Bourdin & Idczak (2014), we define 
the left- and right-sided Riemann-Liouville fractional integrals of 
order cr > 0 of a function ^ 6 L* as 


I=— r 

r(a) X 


Axq{x){t-xf 


'=-r 

r(a) J, 


d.xq(x){x-tr 


(Bl) 

(B2) 


In the remainder of this section, we will use the result illustrated 
by Bourdin & Idczak (2014) in their Proposition 2 (for a proof see 
Samko et al. 1993) 


r dt(c^,)(o<? 2 (o= r dt9,(o(/p?2)(o. 

a *Ja 


(B3) 


Bl Density 

When considering the isotropic limit of the DF, the integral to be 
solved to calculate the density is: 


T'’= f - k). 

W Jo 


(B4) 


We can use fractional calculus to solve it, by changing variable of 
integration (using x = ^ — k) and by considering that 


qi{x) = Ey(g,x), 

qiix) = 1, 

I'o+qdx) ^Ey[g+fx), 


■\fn 


thus obtaining 


IP = 


Jo 


dxE/g+ fx) = Ey(g+ |,0), 


(B5) 

(B6) 

(B7) 

(B8) 


(B9) 


which was solved by using equation (DIO). 

The density of the anisotropic models is calculated by solving 
the integral introduced in equation (A8). To do this, we can use 
the result shown in equation (B3) by noticing (see equations Dll 
and D13) that 


qi(k) = exp(-kp^). 

(BIO) 

1 

II 

(BID 

,1/2 2£(VIp) 

Ifqdk)^ , 

2/np 

(B12) 

II 

+ 

1 

(B13) 

and by rewriting the integral of equation (A8) as 


XP = J' dk exp {-kp^) (s + f<l>-k)- 

(B14) 


This integral can be solved by parts, by using the expression of the 
derivative of the lower incomplete gamma function (equation D4) 
and by using equation (D17) to obtain 


IP = 


1+p" {i + p^)ng+\) 


(B15) 


The last step can be rewritten as equation (Al 1) by using the recur¬ 
rence relations shown in equations (D8) and (D19). 


B2 Velocity dispersion 


In the isotropic limit of the DF, the velocity dispersion is calculated 
by means of an integral with the same structure as the one found 

in equation (B4). The velocity dispersion of the anisotropic mod- 

2 

els is calculated by solving the integral 1'’°' introduced in equa¬ 
tion (A13). We can use the result shown in equation (B3) also in this 
case, by considering the function q^ introduced in equation (BIO), 
its fractional integral (equation B 12), and 

q2{k) = kEy(^g,J>-kj, (B16) 

=^Ey{g+lJ,-k) + kEy{g+f3>-k), (B17) 

? 

and by rewriting the integral as 


XP"^ - J die exp {-kp^) Ey(^g + |, 0 - k) 

-I-2 J dkex^[--kp^^kEy[g + fej)-k^. 


(B18) 


The first integral is in the same form as the one we found for the 
density, and the second one can be reduced to something similar 
with an integration by parts. By solving the integrals, we obtain: 



20^+iiF,(l,g-t |, 




r(g+i) 


g + p <!> + 


1 -t- 


(B19) 


which can be rewritten as equation (A 16) by using the recurrence 
relations shown in equations (D8), (D18), and (D19). 

By inspecting equations (A17) and (A20), it is immediate to 
notice that the radial and tangential components of the velocity dis¬ 
persion are calculated with integrals that can be written as a com¬ 
bination of those solved in this section by means of fractional cal¬ 
culus. 


APPENDIX C: DIFFERENTIAL ENERGY DISTRIBUTION 


The differential energy distribution gives the amount of mass per 
units of energy (Binney & Tremaine 1987). Here we briefly recall 
how to calculate it for isotropic and anisotropic systems. 

For a DF that only depends on E, the differential energy dis¬ 
tribution can be calculated as: 


dM 


(Cl) 


where f{E) is the DF, and g(E) is the density of states, which is the 
volume of phase space per unit energy and is defined as 




g(£)= I d^rdJv6(E-H\ 


(Cl) 


where S(x) is the Dirac delta function. For a spherically symmetric 
systems, this integral can be expressed as 


g(£’)=16;r^ J' drd J' dvv^6^- 


lv^ + (f>-E 


(C3) 


where r^(E) is the radius at which 0 = £. By changing variable 
(using y = 12), we finally obtain 

(C4) 


dr A dl(E - 0), 

0 
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and the differential energy distribution is therefore 

— ^l6n-f(E)j^ drd - cj,). (C5) 

When considering anisotropic systems, for which the DF de¬ 
pends also on the angular momentum J, the differential energy dis¬ 
tribution is obtained as 

dM r , , 

— = \ d\d\6(,E-H)f{H,J). (C6) 

This integral can be expressed as 


^ = f J dVr dv, v, <5(E - J), (Cl) 

and it can then be rearranged by changing variable and introducing 
y = r Vt in the following way 


dM 

d£ 


?,71‘ 




AJ J 6(E - H)f(HJ). 


(C8) 


The integral over is solved by using the fact that 


72 

: 2(£ - O) - 


(C9) 


to obtain 


dM 

d£ 


16 k' 




Jf(E,J) 


d2(E - (F) - 


(CIO) 


By using the expression for the DF of the models presented in this 
paper (equation 1) we can perform the integration over J in this last 
equation and write the differential energy distribution as a function 
of the part of the DF that depends on energy only, f(E): 


dM 

d£ 


16k- f(E) 


J rrm(E) 

0 ‘ 


dr V 2 r^srF 


'dE — (p 


(Cll) 


where F(x) is the Dawson integral (see Appendix D3). In the limit 
of Ta —» CO this reduces to the result for the isotropic case shown in 
equation (C5), which follows from substituting the leading term of 
equation (D14) in equation (Cll). 


APPENDIX D: USEFUL PROPERTIES OF 
MATHEMATICAL FUNCTIONS 

DI Useful properties of the gamma functions 

The gamma function of a positive integer n is defined as 

F(n) = («-!)!, (Dl) 


while for non-inleger arguments a, it can be written as an integral 

(D2) 


F(a) = r Atf * exp(-t). 

Jo 

The lower incomplete gamma function is given by 


y(a, x) = j dlf ' exp(-f), 
Jo 

and its derivative is 
dyffl, x) 


dx 


= x“ ' exp(-x). 


(D3) 


(D4) 


D2 Useful properties of the Ey(a, x) function 

The exponential function Ey(a, x) is defined as 

1 G 

Ey(a,x)= - I drt““* exp(x - r), (D5) 

r(a) Jo 

and an alternative expression is given by means of the lower incom¬ 
plete gamma function 


Ey(a, x) = 


exp(x) 7 (a, x) 


The series representation of this function is 


y(a, x) = y - 

^ Zj Hi -H a 


F(i -t- fl -t- 1) 

The following recurrence relation holds 
Ey(a, x) = Ey(a + l,x)+ — --. 

F(fl -I- 1) 

The derivative and the integral of Ey(a, x) are given by 
dEy(a, x) 


f 


, =Ey(a-l,x), 

dx 

dx Ey(a, x) = Ey(a + 1, x) -t constant. 


(D6) 

(D7) 

(D8) 

(D9) 

(DIO) 


A proof of these equations can be easily obtained by writing 
Ey(a, x) as in equation (D6), and by considering equation (D4) and 
the recurrence relation (equation D8). The convolution formula 


-I— r dyEyia,x-y)y'’'=Ey(a + b,x) 
r(i) Jo 


(Dll) 


can be obtained by using the series representation of Ey(a, x) (equa¬ 
tion D7) and by changing variable, to express the integral with a 
form that allows us to recognize the Beta function: 


B(m,n) = r' dy(l -y)'"-‘y"-‘ = 

Jo r(m -I- n) 


(D12) 


The identity of equation (Dll) accounts for the fractional integra¬ 
tion of Ey(a, x) (see equation Bl). 


D3 Useful properties of the Dawson integral 

The Dawson integral (sometimes called the Dawson function) is 
defined as 


F(x) = exp(-x^) I dyexp(y^). 
Jo 


(D13) 


It is also possible to express F(x) as a sum as 
” (-1)'x2'+>f(|) 

The Dawson integral is an odd function, and its derivative is 
dF(x) 


dx 


= 1 - 2xF(x). 


(D14) 


(D15) 


D4 Useful properties of the confluent hypergeometric 
function 

The confluent hypergeometric function is defined as 
F(a -I- i) r(b) x' 


iFi(a,b,x) =^- 


F(fl) r(b + i) F(i -I- 1) 


(D16) 
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It can also be defined by means of an integral, as 

iFi(a,b,x) = -- r dy exp(j:y)/^‘(l(D17) 

r(a)r(Z? - a) Jo 

which holds for Re(fo) > Re(a) > 0. The following recurrence rela¬ 
tions hold 


iFi(2,b,x) = (2-b + x)iFi(l,b,x) + b-l, 
xtFi{l,b+ l,x) = btFi{l,b,x)-b. 


(D18) 

(D19) 


We also note that this function is related to the exponential func¬ 
tion, and for b = av/e have ]Fi(a,a,x) = exp(x). Another useful 
property of this function is that 

iFi(a,b,0) ^ 1. (D20) 


The density (equation 11) and the velocity moments (equa¬ 
tions 15 - 17) of anisotropic models are expressed by means of the 
function iFi(a, b, x) with a = 1. When considering a = 1 in equa¬ 
tion (D17), we obtain 

iFi{l,b,x) = exp(x) 7 (Z?- l,x). (D21) 

We point out that when x < 0 two of the quantities appearing in 
equation (D21) are imaginary. This is the reason why in general 
this expression cannot be used to speed up the code by express¬ 
ing the equations mentioned above in a more compact way. When 
considering integer values of b, however, we can simplify the hy¬ 
pergeometric function, and express it by means of exponentials and 
polynomials. The smallest value of b we use is g -l- |, therefore b 
assumes the smallest integer value when g = ^ for = 3 we calcu¬ 
late 

,Fi(l,3,x)= 4[exp(x)-l-x]. (D22) 

Combined with the recurrence relation mentioned earlier, the re¬ 
sults of the anisotropic models and half-integer values of g can be 
expressed by means of these elementary functions. 

The asymptotic series expansion for the confluent hypergeo¬ 
metric function when |x| —» oo is given by 


r(^) 

l + O 


Fi(a,b,x) oc ^^( x) 


+ ^ exp(x)x“-'’ 

l + O 

(‘)l 

r(fl) 


U/J 


(D23) 


This is useful to compute the density and velocity moments, be¬ 
cause for very large |x| the evaluation of this function becomes 
inaccurate in python (see section 4.1). In particular, the func¬ 
tions that are needed to compute these quantities, iFJl, 7>, -x) and 
1 ^ 1 ( 2 , b, -x), have the following behaviour for large |x|; 


lim iFJl, b, -x) 


lim 1 ^ 1 ( 2 , b, -X) 


b-\ 

X 

{b-\)(b-2) 


(D24) 

(D25) 
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